A couple of weeks ago, I published the second entry of the Replacing MATLAB with Python series I started in August. I am sharing examples of how to do that, with some comments about the problems I find while learning Python.

In the previous post of the series, I didn’t really share any Python code, though. It was an example of how I would use MATLAB to get information from raw data collected from a race car’s CAN bus.

In this entry,  I’m using Python to do the first part of the task: visualize the signals we have by creating a figure with subplots, filter the raw acceleration data, compute the vehicle speed for the two drivers, and compare it. It took me more than twice the time to do it in Python than it did with MATLAB. That’s mainly because I have been using MATLAB for many years now, and I am just getting used to Python. Especially when formatting plots, there are several things to handle differently from Matlab,  or even from one Python library to another.

As always, bear in mind that this is not intended to be a Python course, but rather a series of practical examples.

Alright, to the code!

Importing the data

I’m following the same workflow that I used for the MATLAB version of this, where I also provided more context on what we have. Short story, we want to compare the performance of two drivers, and we have the car data for each one in separate CSV files.

RaceCar data - Google Drive

Google Drive

This is the code I wrote to start with. We’ve already seen how to load information from CSV files and use Pandas to handle our data.

When using Python, the first thing we do is import the libraries we are going to use. We don’t necessarily have to put the imports on top of everything, but it’s in general good practice.

#%% Import libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

#%% Load data from CSV files

df1 = pd.read_csv('driver1.csv')
print(df1.columns)

df2 = pd.read_csv('driver2.csv')
print(df2.columns)

#%% Define constants
gear_ratio = 16     # Transmission gear ratio
r = 0.263           # Wheel radius [m]

>> Index(['Time_s', 'BrakeFront', 'AccX', 'spd_rpmFL', 'spd_rpmFR', 'spd_rpmRL',
       'spd_rpmRR'],
        dtype='object')

All good! We have a unique time base and 6 signals in each file. But we have to make sure they are consistent and can be used together. This is why the next thing I did in MATLAB was plotting the variables.

When I tried to do the same thing using Python, it took me a while. More than I am comfortable admitting. I was stubbornly trying to get the exact same look I got with the default settings in MATLAB and, obviously, that required a lot of formatting. At some point I got it, but the number of trials with format options was just not worth it. Just look at the plt.rcParams.update call at the bottom of the next code snippet! My advice: if you are replacing MATLAB with Python, don’t try to get exactly the same look on your charts unless there is a real need for that- in my case, there wasn’t.

# If needed, reset plot parameters from any previous plot
# plt.rcParams.update(plt.rcParamsDefault)

fig = plt.figure(figsize = (1200/300, 600/300), dpi = 300)

c = 1 # subplot counter
for col in df1.columns[1:]:
    ax = plt.subplot(3, 3, c)
    ax.plot(df1['Time_s'], df1[col], label = '1', linewidth = 0.2)
    ax.plot(df2['Time_s'], df2[col], label = '2', linewidth = 0.2)
    leg = ax.legend(edgecolor = 'black',
                    borderpad=.5, borderaxespad = .5,
                    fontsize = 3, fancybox = False,
                    loc = 'upper right')
    plt.xlabel('Time [s]', fontsize = 3.5)
    plt.title( col, fontweight="bold", y = 1.05, fontsize = 3.5)
    leg.get_frame().set_linewidth(.1)    
    c = c + 1

plt.rcParams.update({'axes.grid': False,
                    'xtick.labelsize': 3,
                    'ytick.labelsize': 3,
                    'axes.labelpad': 1.0,
                    'axes.labelsize': 3.5,
                    'axes.linewidth': 0.5,
                    'ytick.major.pad': 1,
                    'ytick.minor.left': False,
                    'xtick.minor.bottom': False,
                    'ytick.major.size': 1,
                    'ytick.major.width': .2,
                    'xtick.major.pad': 1,
                    'xtick.major.size': 1,
                    'xtick.major.width': .2})

fig.tight_layout()

fig.subplots_adjust(
    top=0.92,
    bottom=0.1,
    left=0.05,
    right=0.98,
    hspace=0.6,
    wspace=0.2
)

fig

Some of the useful resources I found online to learn about formatting are: first of all, Matplotlib's official documentation. Second, the countless number of sites with Python tutorials online - just googled my questions and found good answers in the first page of results. Third, previously-answered questions on stackoverflow.

Something nice about Python is that we are able to iterate directly through the columns of our dataframe, instead of creating an index and then slicing the dataframe by column position. By doing it this way, I just needed to make a counter to move to the next plot, and the data was selected from the current column. The results are not identical to the ones I got in MATLAB, but that’s OK.

All this additional formatting is needed only if you are using the relatively low-level library Pyplot and you want to get some really specific format. If you loosen the format a bit and open up to other options - like Seaborn for example - you can actually get away with very little code to get some impressive plots. Or, if you need to make interactive plots, you could use Plotly, like I did when visualizing how warped the bed of my 3D printer was. I’ll show another example next, when we compare the speeds of the two drivers.

Fine, we can make plots without stressing too much about the format and live with it.

Let’s move on.

Filtering acceleration data

As in the MATLAB example, I used a digital filter to smooth the acceleration signals. I’m filtering only those because they are the only ones containing raw data from sensors. Here’s what the raw acceleration data looks like when compared against the filtered one:

The code I used for doing this is quite simple:

# Correct and filter longitudinal acceleration data

# Compute the sum of the in-wheel motor speeds
df1['SumSpeeds'] = df1.iloc[:,3:].sum(axis=1)
df2['SumSpeeds'] = df2.iloc[:,3:].sum(axis=1)

# Remove the DC offset from the raw acceleration signals
df1['AccX_NoOffs'] = df1['AccX'] - df1[df1['SumSpeeds'] == 0]['AccX'].median()
df2['AccX_NoOffs'] = df2['AccX'] - df2[df2['SumSpeeds'] == 0]['AccX'].median()

# Filter the corrected signals
from scipy.signal import savgol_filter

df1.AccX = savgol_filter(df1.AccX_NoOffs, 101, 3)
df2.AccX = savgol_filter(df2.AccX_NoOffs, 101, 3)

# reset figure parameters if previously changed
#plt.rcParams.update(plt.rcParamsDefault) 

# Compare raw acceleration vs filtered
fig = plt.figure(figsize = (1200/300, 400/300), dpi=300)
plt.plot(df1.Time_s, df1.AccX_NoOffs, label = 'Raw (removed DC offset)', linewidth = .5)
plt.plot(df1.Time_s, df1.AccX, label = 'Filtered', linewidth = 1)
plt.xlabel('Time [s]')
plt.ylabel("Acceleration [G-units]")
plt.xlim(600, 650)

leg = plt.legend(fancybox = False, 
                edgecolor = 'inherit', fontsize = 4,
                loc = 'upper right')

plt.title('Raw vs Filtered Longitudinal Acceleration')

plt.rcParams.update({'font.size': 5,
                    'axes.grid': True,
                    'axes.grid.which': 'both',
                    'xtick.minor.visible': True,
                    'ytick.minor.visible': True,
                    'grid.linewidth': 0.2,
                    'lines.linewidth': 0.2})

# Display figure
plt.tight_layout()
fig

Implementing the same type of filter as in MATLAB just requires an additional library import. But it is free :). One thing to note is that the input arguments for the Python version savgol_filter(data_array, window_size, poly_order) are swapped with respect to MATLAB’s (recall I used sgolayfilt(data_array, poly_order, window_size)  in the previous post).

One thing many people like about MATLAB (including me) is its documentation. It’s very good! If you need to know something, you type either help or doc  if you want more details, and get the information right there. With Python, you can do something similar using print(<library,object or function>.__doc___) or simply help(<library,object or function>) .

Unfortunately, the official documentation for many Python functions is not as detailed as the MATLAB counterpart. For instance, compare the MATLAB and Python (Scypi) official documentation for the filter we used. The former is much clearer. This can be a big thing for many people wanting to learn Python coming from MATLAB, but there is a way to go around it! There are a lot more tutorials online about Python than there are for MATLAB. You just have to embrace the power of googling, and lose the fear of sites like stackoverflow. For example, I found this great article about the filter we are using. Additionally, the best thing to do is check some Signal Processing book to learn more about digital filters.

Ok! I am deviating from the point of this article. Going back to the code examples.

Finding the vehicle speed: indexing dataframes

As I mentioned in the MATLAB example, we can estimate the vehicle speed using the information from the speeds of the four in-wheel motors. We just need to treat differently the cases in which the vehicle was accelerating and those in which it was braking. To do this, we'l need to index (slice) the dataframe based in conditions. See how I did that in the code below:

import timeit

%%timeit

df1['Vehicle Speed'] = 0
df2['Vehicle Speed'] = 0

# Estimated vehicle speeds when braking (AccX < = 0)
df1.loc[df1['AccX'] <= 0, 'Vehicle Speed'] =  df1.loc[df1['AccX'] <= 0,
                                            ('spd_rpmFL','spd_rpmFR','spd_rpmRL','spd_rpmRR')].max(axis=1)*np.pi/30*r/gear_ratio*3.6
df2.loc[df2['AccX'] <= 0, 'Vehicle Speed'] =  df2.loc[df2['AccX'] <= 0,
                                            ('spd_rpmFL','spd_rpmFR','spd_rpmRL','spd_rpmRR')].max(axis=1)*np.pi/30*r/gear_ratio*3.6

# Estimated vehicle speeds when accelerating (AccX > 0)
df1.loc[df1['AccX'] > 0, 'Vehicle Speed'] =  df1.loc[df1['AccX'] > 0,
                                            ('spd_rpmFL','spd_rpmFR','spd_rpmRL','spd_rpmRR')].min(axis=1)*np.pi/30*r/gear_ratio*3.6
df2.loc[df2['AccX'] > 0, 'Vehicle Speed'] =  df2.loc[df2['AccX'] > 0,
                                            ('spd_rpmFL','spd_rpmFR','spd_rpmRL','spd_rpmRR')].min(axis=1)*np.pi/30*r/gear_ratio*3.6
                                            
# Output 62.2 ms ± 3.78 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

That was very quick. As you can see, if you get used to indexing MATLAB tables, the switch to Pandas dataframes is a lot easier. The syntax (or rather, the logic) is very similar - though it is not exactly the same.

Now, we can just compare the speeds. First, I computed the start and end times for both drivers. This is not precise, I just defined a speed threshold of 20 km/h and found the instants when the car was slower and faster than that. Again, dataframe slicing is very helpful for this:

# Find index of first instant at which speed is greater than 20 km/h
index_start1 = df1[df1['Vehicle Speed']>20].index[0]
index_start2 = df2[df2['Vehicle Speed']>20].index[0]
time_start1 = df1['Time_s'][index_start1]
time_start2 = df2['Time_s'][index_start2]

# Find index of last instant at which speed is greater than 20 km/h
index_end1 = df1[df1['Vehicle Speed']>20].index[-1]
index_end2 = df2[df2['Vehicle Speed']>20].index[-1]
time_end1 = df1['Time_s'][index_end1]
time_end2 = df2['Time_s'][index_end2]

# Compute elapsed times
t_elap_1 = time_end1 - time_start1
t_elap_2 = time_end2 - time_start2

Next, I made a comparative plot (subplots again).  This time with a different approach, though

# %% Compare the speeds of the two drivers
from plotly.subplots import make_subplots
import plotly.graph_objects as go


# Create figure with subplots and add titles
fig = make_subplots(
    rows=2,
    cols=1,
    shared_xaxes=True,
    shared_yaxes=False,
    subplot_titles=((f"Driver 1 approximate time : {t_elap_1:.1f} seconds<br>" 
                     f"Max. speed: {df1['Vehicle Speed'].max():.1f} km/h, "
                     f"mean speed: {df1['Vehicle Speed'].mean():.1f} km/h"), 
                    (f"Driver 2 approximate time : {t_elap_2:.1f} seconds<br>"
                    f"Max. Speed: {df2['Vehicle Speed'].max():.1f} km/h, "
                    f"mean speed: {df2['Vehicle Speed'].mean():.1f} km/h")
))

# Plot data for the first driver
fig.add_trace(
    go.Line(y=df1['Vehicle Speed'] ,x=df1.Time_s,name = 'Driver 1'),
    row = 1, col = 1)
# Add vertical lines to mark start and end points
for time in [time_start1, time_end1]:
    fig.add_shape(go.layout.Shape(type = "line",
                                    yref = "y", xref = "x",
                                    x0 = time, y0=0,
                                    x1 = time, y1 = 100),
                                    row=1, col=1)
fig.update_yaxes(title_text="Speed [km/h]", range=[0, 100], row=1, col=1)

# Plot data for the second driver    
fig.add_trace(
    go.Line(y=df2['Vehicle Speed'] ,x=df1.Time_s,name = 'Driver 2'),
    row = 2, col = 1)
# Add vertical lines to mark start and end points
for time in [time_start2, time_end2]:
    fig.add_shape(go.layout.Shape(type = "line",
                                    yref = "y", xref = "x",
                                    x0 = time, y0=0,
                                    x1 = time, y1 = 100),
                                    row=2, col=1)
fig.update_yaxes(title_text="Speed [km/h]", range=[0, 100], row=2, col=1)

fig.update_layout(title = {
                            'text': '<b>Driver time comparison</b>', 
                            'xanchor' : 'center',
                            'x' : 0.5, 'y': .95, 'font_size' : 20
                            },
                showlegend=False)

fig

Notice that I used Plotly instead of Pyplot. Confusing names, right? Pyplot is the “vanilla” one that comes with Matplotlib. It is easier to learn when starting with Python from a Matlab background and is used a lot for scientific and engineering publication-grade figures. But Plotly offers the possibility to make interactive plots with a little coding, which is awesome. Just look at the results:

Wrapping up

We completed the first part of the task, obtaining the same results as the Matlab example - fortunately! otherwise, there would be a great error somewhere.

When writing this, I think I spent 80% of my time figuring out how to make the plots look like I wanted. The rest was not really that hard. This is not surprising, since Python has gained popularity precisely for being easy to learn. It is easy, as long as you use it as it is intended to. The best way to achieve this is to follow a set of guidelines for coding style, which are known as the “Pythonic” way. I will try to learn this and apply it to my code as I go. I think it makes you think more in terms of code efficiency and make better use of the language.

I hope these code examples were useful for you. In the next post, I will complete the second part of the task, using histograms and statistics to analyze the vehicle speeds and the braking performance of the two drivers.

In the meantime, I hope you have a wonderful day!

It’s been a couple of weeks since I published the initial post of a series about replacing MATLAB with Python. If you haven’t read it yet, the short story is that I forced myself to use Python to do some of the tasks I would normally do using MATLAB at work. I committed to sharing what I've learned in the process, including simple code examples to illustrate it, which may be useful if you are interested in doing something similar. I'm no expert, though, so don't expect to see super complex codes or analyses here!

However, I figured that it would be easier to share separately the MATLAB examples and the Python ones. Mainly because, otherwise, the articles would be excessively long.

In this second entry, I am sharing an example of the workflow I would follow for getting insights from real-world data using histograms and some very basic statistics.

We will be using two datasets containing raw data logged from the CAN bus of a Formula Student Electric car. The data is already stored in a .CSV format to make this more “universal”, since there are different formats for storing automotive data.

RaceCar data - Google Drive

Google Drive

For context, the datasets were recorded while two drivers tested exactly the same vehicle in exactly the same circuit. On the same day, a few minutes apart from each other. There was no transponder installed for recording the lap time accurately, and the information from the electric motors is not available for some reason. This is about all the context we need. We want to analyze the data from both drivers but won’t go into the details.

Let’s jump right into it!

When using MATLAB, I load the data from the CSV files, store the data in some variables, and define the constants that I am going to use later.

%% Get descriptive statistics from CAN datalogs
close all, clc

filename1 = 'driver1.csv';
filename2 = 'driver2.csv';

% Read file 1
t1 = readtable(filename1);
% Verify that I get all the columns I need
t1.Properties.VariableNames

% Repeat for file 2
t2 = readtable(filename2);
t2.Properties.VariableNames

%% Define vehicle parameters
% This was a Formula Student Electric car, with four in-wheel electric
% motors coupled to an epicyclic transmission, all four with the same gear ratio
gear_ratio = 16;    % Transmission gear ratio
r = 0.263;          % Wheel radius [m]

I'm using tables to handle the data because I find them very easy to use. Additionally, using tables in MATLAB can help you transition more easily to using Pandas data frames if you want to. This will be more evident later on. Tables are not the most efficient way to handle very large datasets, but they do offer some convenience and if you know how to index them, they work great.

Let’s move on. I have the data loaded on MATLAB, so now I can do stuff. As I’ve mentioned before, a good practice is to start with a simple plot or, in this case, a bunch of them,  since we have few variables and files.

%% Plot the data available, just to have a quick overview
figure('Position', [10 10 1200 400])
tiledlayout('flow', 'TileSpacing', 'compact')

for i_plot = 2:size(t1,2)
    
    nexttile
    plot(t1.Time_s, t1{:, i_plot}, 'DisplayName', '1')
    hold on
    plot(t2.Time_s, t2{:, i_plot}, 'DisplayName', '2')
    title(t1.Properties.VariableNames{i_plot}, 'Interpreter' , 'none')
    legend
    
end

Great! we have comparable data in both files. We have pressure data from a sensor in the hydraulic circuit of the car, raw acceleration data from an IMU, and the rotational speed of the four motors attached to the wheels. There is a slight time offset between the two data sets, and one of them is a bit longer. This is raw data, so it is something completely normal.

Let’s find out which driver was faster.  To do this, we need to compute the vehicle speed.

Estimating the vehicle speeds

This is not easy at all - it might seem, but it is not - at least if you want to do it accurately. What we can do is settle for a simplified calculation assuming the following:

• When the vehicle is accelerating, the speed of the car is similar to the speed of the contact point of the slowest tire, i.e. the one that rotates slower than the others and thus has a smaller slip ratio.
• On the other hand, when the car is braking, the slowest tire is most likely locked, with a higher absolute slip ratio than the others. So in this case the vehicle's speed is closer to that of the fastest spinning tire - no accelerating torque is present.

This is a rough estimate but should be enough for this analysis. But before we do this, we need to do some work to clean up the raw acceleration data. I will remove the DC offset, and filter the noise.

%% Filter the longitudinal acceleration

% Compute the sum of the wheel speeds
sum_of_speeds1 = t1.spd_rpmFR + t1.spd_rpmFL + t1.spd_rpmRR + t1.spd_rpmRL;
sum_of_speeds2 = t2.spd_rpmFR + t2.spd_rpmFL + t2.spd_rpmRR + t2.spd_rpmRL;

% Remove any DC - offset from the raw acceleration data
t1.AccXNoOffs = t1.AccX - median(table2array(t1(sum_of_speeds1 == 0, 'AccX')));
t2.AccXNoOffs = t2.AccX - median(table2array(t2(sum_of_speeds2 == 0, 'AccX')));

% Filter the signals
t1.AccX = sgolayfilt(t1.AccXNoOffs, 3, 101);
t2.AccX = sgolayfilt(t2.AccXNoOffs, 3, 101);

% COmpare the signals
figure('Position', [10 10 1200 400]);
hold on
plot(t1.Time_s, t1.AccXNoOffs, '-k', 'DisplayName', 'Raw AccX (removed DC offset)', 'linewidth', .5)
plot(t1.Time_s, t1.AccX, '-r', 'DisplayName', 'AccX (filtered)', 'linewidth', 1.5 )
xlim([600, 650])
xlabel('Time [seconds]')
ylabel("Acceleration [G - units]")
legend
grid minor

I filtered the acceleration signals using a Finite-Impulse-Response filter, in this case, a Savitzky–Golay filter. This type of filter does a good job preserving the height and location of the peaks and the area below them. Implementing this type of filter in MATLAB is as simple as just typing the sgolay command and including some parameters chosen accordingly. But only if you have the Signal Processing Toolbox license. So this might be one of the reasons to consider using free software options. Of course, with a little more work, you can implement a similar filter yourself using “standard” MATLAB commands. (edit: if you actually check the book Introduction to Signal Processing cited in the MATLAB documentation of the filter, you can actually see the implementation for many types of digital filters both in MATLAB and C, and download the files. Cool!)

Now, how to determine the vehicle speed using the approach I mentioned before?

Easy, right? Just make a for loop and iterate through the rows of the table. When the acceleration is positive take the minimum wheel speed, when it is negative, take the maximum one. And of course, convert each value to linear speed.

%% Estimate the vehicle speed: the for-loop way
% Extract the wheel rotational speed data from the two files
wRPM1 = t1(:,{'spd_rpmFL', 'spd_rpmFR', 'spd_rpmRL', 'spd_rpmRR'}); 
wRPM2 = t2(:,{'spd_rpmFL', 'spd_rpmFR', 'spd_rpmRL', 'spd_rpmRR'});

tic
for i_ax1=1:height(wRPM1)
    
   if t1.AccX(i_ax1,1)>0
       VEL1_for(i_ax1,1)=min(wRPM1{i_ax1,:})*(pi/30)*r/gear_ratio*3.6;
   else
       VEL1_for(i_ax1,1)=max(wRPM1{i_ax1,:})*(pi/30)*r/gear_ratio*3.6;
   end
end
toc % Elapsed time is 289.592689 seconds. <--- NOT GOOD!

That took forever! I don't want to wait that much time when doing this kind of thing! What if the data grows larger, or if I need to compare 100 files instead of 2, should I wait 4 hours? Maybe my PC is not the fastest, but there has to be a better way!

Fortunately, yes, there is a better way. We can use indexing methods instead of for loops.

%% Estimate the vehicle speed: Release the power of table indexing
tic

% Create empty columns for storing the vehicle speed
t1.Veh_Spd = zeros(height(t1), 1);
t2.Veh_Spd = zeros(height(t2), 1);

t1(t1.AccX > 0, 'Veh_Spd') = table(min(t1{t1.AccX > 0,...
    {'spd_rpmFL','spd_rpmFR','spd_rpmRL','spd_rpmRR'}}, [], 2)*pi/30*r/gear_ratio*3.6);
t1(t1.AccX <= 0, 'Veh_Spd') = table(max(t1{t1.AccX <= 0,...
    {'spd_rpmFL','spd_rpmFR','spd_rpmRL','spd_rpmRR'}}, [], 2)*pi/30*r/gear_ratio*3.6);     

t2(t2.AccX > 0, 'Veh_Spd') = table(min(t2{t2.AccX > 0,...
    {'spd_rpmFL','spd_rpmFR','spd_rpmRL','spd_rpmRR'}}, [], 2)*pi/30*r/gear_ratio*3.6);
t2(t2.AccX <= 0, 'Veh_Spd') = table(max(t2{t2.AccX <= 0,...
    {'spd_rpmFL','spd_rpmFR','spd_rpmRL','spd_rpmRR'}}, [], 2)*pi/30*r/gear_ratio*3.6);
    
toc % Elapsed time is 0.018584 seconds.

Awesome! we get both velocity vectors in just a fraction of a second, instead of only one in 5 minutes. I just took advantage of the indexing capabilities of MATLAB - it was designed to work with matrix operations, so why not use them?  Let me explain a bit:

I first added empty columns to the tables. Then, I indexed the tables, selecting only the rows for which the acceleration was positive, and the empty column I just created.

In that selection, I inserted the minimum of the motor speeds for each row, again, considering only the rows for which the acceleration was positive, and converting that to linear speeds.

Then I did the same for the opposite case when the acceleration was negative.

A similar approach can be used to avoid for loops when using numeric arrays or cell arrays instead of tables. The difference is just the syntax for indexing them.

Moving on...

Comparing the speeds of both drivers

Let’s see what we’ve got as vehicle speeds:

%% Create a figure with tiled layout
                                                
time_series = figure('Name', ['Time series comparison'], 'Position', [10 10 1200 400]);
tiles = tiledlayout(2,1);
title(tiles, 'Driver time comparison')
xlabel(tiles, 'Time [s]')

%% Plot data from the first file

nexttile
plot(t1.Time_s, t1.Veh_Spd, '-b')
xlim([0, 1200])
ylim([0, 100])
ylabel('Vehicle speed [km/h]')
hold on
time_start1 = t1.Time_s(min(find( t1.Veh_Spd > 20)));
time_end1 = t1.Time_s(max(find(t1.Veh_Spd > 20)));
t_elap_1 = time_end1 - time_start1;

plot([[time_start1 time_end1]; [time_start1 time_end1] ], ...
    repmat(ylim',1, 2), 'Color', 'k', 'LineWidth', 2);
title(['Driver 1 approximate time: ' num2str(time_end1 - time_start1, '%0.1f') ' seconds'...
       newline 'Max. Speed: ' num2str(max(t1.Veh_Spd), '%0.0f') 'km/h, '...
       'Mean Speed: ' num2str(mean(t1.Veh_Spd, 'omitnan'), '%0.0f') 'km/h'])
grid minor


%% Plot data from the second file

nexttile
plot(t2.Time_s, t2.Veh_Spd, '-b')
xlim([0, 1200])
ylim([0, 100])
ylabel('Vehicle speed [km/h]')
hold on
time_start2 = t2.Time_s(min(find(t2.Veh_Spd > 20)));
time_end2 = t2.Time_s(max(find(t2.Veh_Spd > 20)));
t_elap_1 = time_end1 - time_start1;
        

plot([[time_start2 time_end2]; [time_start2 time_end2] ], repmat(ylim',1, 2), 'Color', 'k', 'LineWidth', 2);
grid minor
title(['Driver 2 approximate time: ' num2str(time_end2 - time_start2, '%0.1f') ' seconds'...
       newline 'Max. Speed: ' num2str(max(t2.Veh_Spd), '%0.0f') 'km/h, '...
       'Mean Speed: ' num2str(mean(t2.Veh_Spd, 'omitnan'), '%0.0f') 'km/h'])

Ok, the second driver was a lot faster. There could be very noticeable differences in the driving styles since they were using the same car. So let’s see If we can spot them. First, let’s look at their speeds from a statistical perspective:

%% Get statistics about the driving speeds                                               
histograms = figure('Name', ['Histograms comparison'], 'Position', [10 10 1200 400]);
tiles = tiledlayout(2,1);
title(tiles, 'Comparing vehicle speed data using histograms')
xlabel(tiles, 'Vehicle speed [km/h]')
binEdges = 0:5:100;

nexttile
VEL1 = t1.Veh_Spd;
VEL2 = t2.Veh_Spd;
n_vel=max(numel(VEL1),numel(VEL2));

VEL1(end+1:n_vel)=nan;
VEL2(end+1:n_vel)=nan;
% "Old" method for creating a histogram
hist([VEL1, VEL2], binEdges);
ylabel('Counts')

legend('Driver 1', 'Driver 2')
title('Bar plots of histogram data (hist command)')
grid on

nexttile
% Newer way of creating a histogram
h1 = histogram(t1.Veh_Spd, binEdges, 'normalization' , 'probability')
hold on
h2 = histogram(t2.Veh_Spd, binEdges, 'normalization' , 'probability')
legend('Driver 1', 'Driver 2')
title('Plotting two overlapping histograms (histogram command)')
ylabel('Probability')
grid on

%% Get some descriptive parameters
v1_mean = mean(t1.Veh_Spd, 'omitnan')       % Means of the two samples
v2_mean = mean(t2.Veh_Spd, 'omitnan')
v1_median = median(t1.Veh_Spd, 'omitnan')	% Medians of the two samples
v2_median = median(t2.Veh_Spd, 'omitnan')
v1_std = std(t1.Veh_Spd, 'omitnan')         % Standard deviations of the two samples
v2_std = std(t2.Veh_Spd, 'omitnan')

The histograms on the right side are normalized so that we can compare both files, even if the time of the second driver was shorter - which means that we have fewer data points. By dividing the counts of each bin by the size of the total sample, we actually get a representation of the discrete probability distribution. In this case, the width of the bins is constant, but if there were bins of different sizes, we would also need to divide by the width of the bins to get a valid comparison.

Analyzing the braking pressures

The histograms show that the second driver spent more time at speeds above 60 km/h than the first one.  But both drivers were using exactly the same car. And the mean and median speed values are not that different between the two. So what could have changed? Perhaps it is also important to see how they were using the brakes.

%% Create histograms of the braking pressures
% This considers only the use of mechanical brakes.
% The electric car can also use the electric motors for regenerative braking. 
% This is enabled through a small portion of the brake pedal travel.
% These two particular datalogs do not contain information about that, so we will
% focus on the mechanical (hydraulic) brakes.

t1.BrakeFront(t1.BrakeFront <=0) = 0; %filter negative pressure values
t2.BrakeFront(t2.BrakeFront <=0) = 0; %filter negative pressure values
%%
figure;
hold on
histogram(t1.BrakeFront, 'Normalization', 'probability')
histogram(t2.BrakeFront, 'Normalization', 'probability')
ylabel('Probability')
xlabel('Pressure [bar]' )
title('Brakes pressure in the front circuit ')
legend('Driver 1', 'Dirver 2')
grid on

The high bar on the left side represents the low-pressure values, which are the most common values present during driving. This makes sense since one doesn’t drive with the brakes pressed, especially during a race. But driver 1 did spend more time braking. We can see that the probability of having the brakes completely disengaged (e.g.about zero hydraulic pressure in the circuit) is lower for driver one. And we also see that the probability of having low pressures in the circuit is higher in case 1. This could mean that the first driver actually had the brakes slightly pressed constantly, perhaps without noticing.

Let’s see also at what speeds the drivers were braking.

Finding the vehicle speeds at the start of the braking maneuvers

To do this, we need to identify the speeds at the start of the braking maneuvers. I defined a small pressure value to use as a threshold for identifying when the brakes are pressed. This is a bit of an iterative process - you don’t get it right the first time and need to adjust the value until the results make sense.

Let's use again the power of indexing, I don't want to wait forever with a for loop. We just need to determine in which rows there was a change in pressure that we could recognize as a pedal press. We also need to consider only the cases in which the vehicle was actually moving. So we have three conditions:

1. Low pressure on instant t
2. High pressure (greater than a threshold) on instant t+1
3. Vehicle speed greater than a minimum value at instant t+1

I indexed the tables using these three conditions and used parentheses to make them more readable (they are not really needed for the code to work). The trick is to impose the second and third conditions (t+1) by indexing the table from the second row. This creates the shift in time values that we would normally obtain using a loop with an increasing counter:

%% Find the index of the start of braking

%Define a pressure threshold value to consider as start of braking
p_start = 2; % bar

% Create indices for the start of the braking maneuvers
id_start1 = (t1.BrakeFront <= p_start) & ([t1.BrakeFront(2:end);0] > p_start) & ([t1.Veh_Spd(2:end); 0] >= 15);
id_start2 = (t2.BrakeFront <= p_start) & ([t2.BrakeFront(2:end);0] > p_start) & ([t2.Veh_Spd(2:end); 0] >= 15);

%% Verify that the identified points make sense
figure
plot(t1.Time_s, t1.BrakeFront, 'b', 'LineWidth', 1,'DisplayName', 'Pressure')
hold on
x = [t1.Time_s(id_start1), t1.Time_s(id_start1)];

% Show a black vertical line at the start of each braking maneuver
line(x, ylim, 'Color', 'k', 'LineWidth', .5, 'HandleVisibility', 'off'); 
yline(p_start, 'DisplayName', 'Brake start threshold' ,'LineWidth', 2) 
ylabel('Pressure [bar]' )
yyaxis right

plot(t1.Time_s, t1.Veh_Spd, 'r', 'LineWidth', 1, 'DisplayName', 'Vehicle Spd')
set(gca, 'Ycolor', 'r')
ylabel('Vehicle speed [km/h]')
xlabel('Time [s]' ), xlim([750, 800])
legend, grid on

Great! we have identified the starts of braking maneuvers. Now let’s look at the speeds at which the drivers started braking. We just have to create histograms of the vehicle speeds using the indices we just found:

%% Create histograms and plot them
                                                
histograms = figure('Name', ['Histograms comparison'], 'Position', [10 10 900 600]);
title('MATLAB figure')
xlabel('Vehicle speed [km/h]')
binEdges = 10:5:100;

histogram(t1.Veh_Spd(id_start1), binEdges, 'Normalization' , 'probability')
hold on
histogram(t2.Veh_Spd(id_start2), binEdges, 'Normalization' , 'probability')
legend('Driver 1', 'Driver 2')
ylabel('Probability')
xlabel('Vehicle speed at the start of braking [km/h]')
grid on

Looking at the histograms, we can see that the first driver initiated more braking maneuvers between 60 and 70 km/h, but he also initiated more braking maneuvers at lower speeds - below 55 km/h. The second driver was braking at higher speeds.

This was a simple toy example of how to use MATLAB to do some tasks with real-world data. We used histograms to get valuable insights from the data, something that would not be easy if we just used time series plots to visualize the two datasets.

Sometimes we can be too concerned with the time-based evolution of the data we are looking at, rather than their distribution, or the trends and correlations between different variables.

Seeing how histograms can be easily used to get insights without actually needing to look at the times series can sometimes open a new perspective when looking at data recorded from physical things (cars, machines, sensors, etc).

Hopefully, this example was useful or interesting for you if you are just getting started doing similar things. I would be glad if it was in any way :)

In the next entry, I will repeat the analysis (at least the main part) using Python instead of MATLAB. Stay tuned!

I work a full-time job as an automotive Control Systems Engineer, and - usually - a big chunk of my time during the day is employed using MATLAB/Simulink for different types of tasks.

I recently spent a couple of weeks without an individual MATLAB license at work due to different issues that I am not going to address. Of course, I had the chance to use workstations with shared licenses in the office to continue working in Model-Based software development. For that kind of activity, you can't just go around the license problem, since you work on shared files and projects within a team, and you need to produce models (code) that are part of a larger project.

However, I took the chance to sharpen my Python skills while completing some of the other tasks I had to do, mainly when analyzing data. I could have used other alternatives, like GNU Octave - which even has the same MATLAB syntax, or Scilab. But I wanted to learn more about Python, since it is a general-purpose programming language.

I figured that I could share some of the things I learned, and started to write an article that suddenly was too long for a single post. So this will be an ongoing series of posts. In each entry, I will share a simple but practical example of how to do a certain task with MATLAB, but using Python instead. Probably some of the code can be improved, but I think it can be useful to see these examples if you are a beginner like me.

I think this is important since Python opens the possibilities for individuals and companies that cannot afford to buy an expensive MATLAB license, or that already using MATLAB but have use cases for which MATLAB can be easily replaced with Python, cutting down the license costs.

MATLAB is still being used widely, especially in the Automotive and Aerospace industries, but the usage of Python has grown dramatically in many different sectors, since it can be used for many other things - like I said, general-purpose. Python cannot really be compared with MATLAB, as they were designed for different purposes, but I will focus on the use cases that can overlap between the two of them.

About Automotive file formats and data extraction

Skip this part if you are not interested in automotive-specific data

Extracting and doing the first-level analysis of data from large log files is a very common use case of MATLAB in automotive testing environments - at least from what I've seen in three different large companies in Europe. It could be seen as a mundane activity, but it is not straightforward if the files you have are in different formats. Also,  it might be something you need to do very frequently.

For this series of posts, I will look at the tasks I learned to do with Python once the data was available in some kind of universal text encoding. For instance a CSV file or an ASCII file. This is not always obvious, and there were many cases in which I needed to convert VECTOR CAN log files or other ASAM standard formats to text, or use some kind of API to be able to work with them in Python or MATLAB. But that will be part of another article later on.

How to start using Python?

First of all: this will not be a Python course. I am assuming you have already some basic understanding and want to see some basic examples of how to use it. You can learn the basics very quickly with tons of free tutorials and courses out there. I followed the Scientific Computing with Python Certification on FreeCodeCamp.org a while ago, and it was great to grab the basics. There are also great resources like TechWithTim.net, and many courses on Coursera, EdX, DataCamp, and others.

Python can be actually written in any raw text editor, like Numpad++ , and then run from the Windows command prompt. But you probably don't want to do that.

If you are coming from MATLAB and want to start using Python, you probably want some kind of user interface or IDE (Integrated Development Environment) to work with. There is no unique or best answer here, so I will just tell you about the ones I have tried.

The first thing is, of course, to download and install Python 3 itself (Python 2 is already in an End-of-Life stage and will no longer be supported by new libraries). Then, you can download and install one of the following IDEs:

1. Jupyter Labs. Probably the most used interface for working with Python nowadays. It allows you to work on your web browser, and the Jupyter Notebooks allow splitting the code into cells that can be run individually, as we do with MATLAB's code sections.
2. SPyder IDE. In my opinion, this one provides the closest experience to MATLAB in terms of interface. You can even have a panel to look at the variables you are working with, for example, and interactive plots. You can install add-ins to work with Jupyter Notebooks within the Spyder IDE too.
3. Visual Studio Code (VS Code) This is a more general-purpose programming IDE, but has great support for Python, including support for Jupyter Notebooks too.

Personally, I prefer to use VS Code since it allows me to work also on other types of programming projects in the same environment, for example when working with Arduino. However, the code examples I will share will work no matter the IDE you are using.

Ok! enough with the introduction, let's jump to the first and most simple example of this series.

Example 1: Simple time series visualization

Normally, working with MATLAB or not, the first thing to do after loading data from some large file would be to plot a time series of the most interesting data (or signals, if we are talking about a CAN data log file, for example). This way I can at least know what I’m dealing with. For this example, I will use this example file, which contains some CAN data from a Formula Student Electric car.

test_FSCar_ddMMYYYY_TrackN_setupID2.csv

Google Docs

If you download the file, you'll notice that there are missing values on each variable (CAN signal). This happens a lot especially if you have CAN data, which travels in packets, sometimes at different rates. When the data is exported to other formats, there are different ways to handle the different rates of the signals. In some cases, you will get an output file containing data with individual time vectors for each signal. In this case, the user preferred to have a single time base in the exported file, common to all signals, but the drawback is that there are empty values. I personally prefer this way, since it is easier to do time-based analyses with a unique time vector. We just need to handle the empty values.

To look at this file with MATLAB, I would do something like this:

Notice that I extracted the indices of the non-empty values of each signal, and then used that to represent the data.  I also did some formatting to the plot to make it understandable.

I needed to represent data containing information about speed, in rpm, and torque, in Nm. And I wanted to visualize both against the same time scale, in order to understand what happened during the test I am looking at.

This is why I created a tiled layout. I could have used subplots instead, and the results would have been very similar. Or I could have plotted both speeds and rpm against a single abscissa, but using two different y axes using yyaxis right.

The output will be something like this:

Simple enough, right? Let's see what we have to do instead when using Python.

The first difference you'll notice is that we need to import the Python libraries we'll use at the start of the code. In this case, we are dealing with data from a text file, so it is a common practice to use the library called Pandas. Although it can be slow for doing complex computations or handling large datasets, there are some workarounds. In general, it is better to do the calculations using another library, Numpy. We simply import those using the import command before using anything from them.

Then we can read the file, and load the data as a Pandas dataframe. In Python, we have to get used to working with methods other than functions. This will allow us to do different things with the variables we are working with, depending on what type of object (class) they are. For instance, we loaded the data from the file as a dataframe object that we called df (we could have called it anything, like Pikachu, if we wanted, but df is more self-explanatory and common practice). Then, to extract the non-empty values we use the . notation to access the notnull() method of the dataframe class.

This is different from the function-based approach that I used in the MATLAB script above, where I used the isnan function, and it would have worked on either a double array, a scalar, a matrix or a multidimensional array.

However, methods and classes are not unique to Python, they exist in almost any programming language, and even in MATLAB. The difference is how they are used depending on the intent of the language. MATLAB is specific for matrix computations, so it relies heavily on functions, whereas Python has specific libraries for different purposes, and they use both functions and methods to accomplish some tasks. There are also some built-in functions in Python

Next, we create the figure and plot the data. We first need to specify the figure size, otherwise, the default figure dimensions will be too small. In this case, we are using subplots to obtain the following result:

Granted, the Python plot needs a bit more manual formatting to have the same look as the MATLAB one. But this is not really an issue once you know what you need to change to get the format you want. You can set the defaults for creating figures on the top of your script, or make a function that does the formatting for new figures.

If you have trouble getting some details to work, you can just google your answer. The amount of Python users is huge and growing, so it is very likely you will find someone already asked a similar question online. For instance, I was having a hard time getting the axes titles to show correctly, and there was a Stackoverflow thread solving exactly the same issue - turns out you need to add plt.tight_layout() for the figure to be padded correctly. Just part of the learning curve. It's like when I needed to call the legend in MATLAB as legend('interpreter', 'none') to avoid tex interpretation of variable names with underscores, for instance.

Plus, this example was done using Pyplot. There are other libraries for plotting, which can produce cleaner plots in different cases, depending on your purposes.

This post was a bit longer than the others, being the introduction of the series, but I hope you found it helpful or at least interesting. Over the next weeks, I will be sharing other example codes, highlighting the differences between MATLAB and Python, and showing also the things I struggled with. I hope you are motivated to learn Python as I am, and this can be an incentive to try it on your own.

I will update this first post to index the new examples I include as I go, so it is easier to navigate between the different entries.

Until the next one!

Luis

In 2019 I finally got my first hobby-grade, FDM 3D printer kit.

It was a Creality Ender 2, a cantilevered version of the successful Ender 3. It was the perfect fit for a small apartment, requiring a very reduced space to be installed on a desktop. I was really excited when it arrived, assembled it in 30 minutes, and started printing right away. The result of the first print turned out to be a decapitated version of the iconic Maneki-Neko cat, since the filament sample that came with the printer was not enough to complete the test .gcode file. I ended up modeling a head to fit the incomplete print and shared it online so that other fellow makers could use it as well.

At the time, all I wanted to do was to learn more about 3D printing, mod and upgrade my small Ender 2 and see how I could get the best out of such a cheap machine. I was still at university, completing my master's degree and I had just started working as a research assistant in the Aerosol Technology Laboratory. I had time to tinker around with the printer, right?

I had no idea what I was getting into.

The printer became a Frankenstein and a learning project for me. I switched the printer's 8-bit mainboard for a 32-bit aftermarket one, upgraded the firmware, and installed a part-cooling fan and motor dampers: My first goal was to make it quieter.

Then I wanted to improve the print quality, so started buying other "budget upgrades",  like a Bondtech BMG knock-off from AliExpress. You see, all the parts I was installing were usually sourced from China, so they were cheap, and took ages to arrive. But buying cheap parts for a cheap printer meant that I wasn't spending much after all... right? After all, it was working just as well as a much more expensive printer. What's a couple of euros here and then to make the machine work better? A new print surface to improve bed adhesion. A new set of PTFE tube couplers because the original ones broke...

But buying cheap parts for a cheap printer meant that I wasn't spending much after all... right?

Reliability

It was all good until I started actually wanting to use the printer more consistently during the 2020 COVID-19 lockdown. I printed a bunch of parts for the laboratory I was working on - since we had the equipment, we rushed to provide pre-certification tests for community face masks. Precisely then, I started having serious problems with the prints not adhering to the print bed. I couldn't launch a print and walk away without having to check every 10 seconds on whether the first layer was OK or not.

I tried every imaginable solution and read every possible Reddit thread or forum entry about similar problems. Tried different filaments, with different temperatures. Tested different combinations of filaments and heated bed surfaces, and re-leveled the bed a million times. Until I figured out that the bed was warped, and I needed to compensate for that. So I installed an Auto Bed-Leveling Sensor (3D touch) which would probe the bed before each print to compensate for the warped bed. It only worked partially. Some sports were too close to the bed, some too far.

On top of that, the extrusion started to be very inconsistent. And for various reasons. Sometimes I would fix things by changing the slicer parameters. But at some point, even slicer changes wouldn't help, and I started to find issues with hardware components. One time it was the hot-end heated cartridge not being able to provide enough power. Later on, the extruder gears started making some noise and I found out they were super worn out.

When I look back at the orders for upgrades and replacement components I made through AliExpress, Amazon, 3D Jake and eBay combined, I realize that the total is € 314! That is almost twice the price I paid for the printer itself, and it's not including the time I spent working on it, and reading online to figure out how to solve the problems I found.

I made a list with all the components i bought for the printer, their price, intended use, and how effective those upgrades were.

Unfortunately, I have lost all the data I collected about noise levels, number of successful first layers, bed level meshes, and so on. But I still think the list can be useful for someone starting with 3D printing as I did.

Ender 2 - Modding and running costs - Google Drive

Google Drive

It all comes down to what your plans are, and how you value your time and money

In retrospect, the process of calibrating and upgrading the old Ender 2 and "fighting" to make it work better was an invaluable learning experience for me. I followed a bunch of useful online guides and learned about how awesome the maker community is online. I now understand a lot more about how 3D printers work: the Marlin software, extruder types, motion systems,  stepper motors, power supply units, and got a lot of experience with slicer parameters. It was really a ton of fun too.

So now, when someone asks me: "should I get a 3D printer? How much should I spend on it?", the answer is, of course, a big: it depends.

What do you want to use it for? If you need consistency out of the box, maybe even a cheap Creality machine will do. If you need reliability and quality out of the box, you'll need a bigger budget for sure.

If you just want the printer for hobby use, then ask yourself: Do I want to tinker with the printer as a hobby, and am I willing to spend time learning about how it works? or do I want to use the printer to produce parts for my actual hobby?

That's where I made the decision to get a new one. Now I just need something to print consistently, since I no longer have all the time I had as a student. Plus, the Ender 2 has been discontinued, which means that replacement parts are no longer being produced. This is especially bad for the warped heated bed which had me struggling for the last couple of months. For reference, look at this plot I made with the data from a G29 command (the glass bed was installed!):

For my current needs, I'd rather have a printer that just works to make stuff for my actual hobbies. I think I've learned the basics quite well now.

By the way, some people were interested in the code I used to make plot shown above, so here it is the python script :)

Hi! As you might have guessed by the site name, my name is Luis and I make stuff when I am not working my full-time office job. I am a Mechanical and Automotive Engineering graduate, but I have always loved making things at home: sometimes just drawings, paintings, or small sculptures. And sometimes I do more technical stuff, like python scripts, small electronic projects with ESP8266 or ESP32, 3D printed objects... you get the idea. I make stuff.

People who have known me at different moments of my life might remember me for different aspects. For instance:

• In school, I used to be the guy doing hand sketches and drawings
• During my first years of university, I was the guy obsessed with 3D CAD modeling
• While completing my master's studies, I was the guy who wanted to learn everything about FEM stress analysis - and actually wrote my thesis almost 80% around that.

But now, 2 years after graduating, I realize that what I really enjoy is learning new stuff and applying it to real problems. If you asked my friends from university almost none of them would think I'd be interested in data analysis, electronics, sensors, IoT, and all that kinds of things. But here I am :)  It's a random mix, for now.

I plan to use this site to share things that I have learned and little projects that I do in my spare time. Hopefully decently documented so you might find something useful or at least interesting. I might also share more engineering-related content about things I've done in the past. There is a double-fold benefit for me: I also plan to use this to organize myself better and have a more structured way to actually complete the projects I start.

Now, from all the things I will post, about 80% are self-learned, and I am in absolutely no way an expert on any of the topics I will discuss here, but at least I will do my best to point you guys out to useful resources I used for learning myself.

I have been meaning to publish this for some time now, and I have finally done it (during a week of COVID-19 home isolation).

Wish me luck! and I hope you have a great day!